A NEW MODEL FOR SHALLOW VISCOELASTIC FLUIDS 

FRANgOIS BOUCHUT AND SEBASTIEN BOYAVAL 

Abstract. We propose a new reduced model for gravity-driven free-surface flows of shallow 
viscoelastic fluids. It is obtained by an asymptotic expansion of the upper-convected Maxwell 
model for viscoelastic fluids. The viscosity is assumed small (of order epsilon, the aspect ratio of 
the thin layer of fluid), but the relaxation time is kept finite. Additionally to the classical layer 
depth and velocity in shallow models, our system describes also the evolution of two components 
of the stress. It has an intrinsic energy equation. The mathematical properties of the model are 
established, an important feature being the non-convexity of the physically relevant energy with 
respect to conservative variables, but the convexity with respect to the physically relevant pseudo- 
conservative variables. Numerical illustrations are given, based on a suitable well-balanced finite- 
volume discretization involving an approximate Riemann solver. 
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1. Introduction: thin layer approximations of non-Newtonian flows 

There are many occurences of free-surface non-Newtonian flows over an inclined topography 
in nature, for instance geophysical flows: mud flows, landslides, debris avalanches .... Their 
mathematical prediction is important, typically for safety reasons in connection with land use 
planning in the case of geophysical flows. But their modelling is still difficult, as one can conclude 
from the continuing intense activity in that area (see the reviews [1] [2D] e.g. plus the numerous 
references cited therein and below). In this paper, our purpose is (i) to derive a new simple model 
for a thin layer of viscoelastic non-Newtonian fluid over a given topography at the bottom when 
the motion is essentially driven by gravity forces and (ii) to numerically investigate the prediction 
of that simple model in benchmark cases. 

Our methodology follows the standard derivation of the Saint- Venant model for gravity-driven 
shallow water flows, as developped in [26J for instance. For non-Newtonian fluids, there already 
exist similar projects in the literature. But to our knowledge, they use different models as starting 
point: power-law and Bingham models [24l [T^, or a kinetic model for microscopic FENE dumb- 
bells HO] (see also the Section 16.21 for comparison with a kinetic interpretation of our model using 
Hookean dumbbells). Here, we derive a reduced form of the Upper-Convected Maxwell (UCM) 
equations, a widely-used differential model for viscoelastic fluids valid in generic geometries, in 
the speciflc case of gravity-driven free-surface thin-layer flows over an inclined topography. In 
particular, the influence of each term in the equations is compared with the aspect ratio 

(1.1) V^«e<l 

between the layer depth h and its longitudinal characteristic length L as a function of a small 
parameter e. 

The reduced model obtained is computationally much less expensive to solve numerically than 
the full UCM model with an unknown free surface (compare for instance with numerical simulations 
in [3T] of a full 3D model). So we can easily investigate its predictions in a number of test cases, 
to show the capabilities of the model. We note two important aspects from the mathematical 
viewpoint to discretize our model. It is endowed with a natural energy law (inherited from the 
UCM model) but has a non-standard hyperbolic structure (the physically relevant energy is not 
convex with respect to the conservative variables). These features of our model have important 
consequences on the numerical simulation. Whereas we can only perform numerical simulations in 
a formal way (because the non-standard hyperbolic structure does not fit in the usual numerical 
analysis), we can nevertheless confirm that they are physically meaningful (owing to the natural 
energy law. satisfied at the discrete level). 

Regarding the literature, we would like to make two further comments in order to better situate 
our work and its originality. On the one hand, numerous models for thin layers of non-Newtonian 
fluids have already been derived in the physics literature. We are aware of only one reduced version 
of the UCM model which is very close to ours, see [2H[13] and a sketch of that work in [32] ■ But the 
reduced model, obtained with another methodology and with a different perspective (ad-hoc model 
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to investigate the break-up and swell of free jets and thin films rather than asymptotic analysis of 
general fluid equations), finally applies in different conditions (without gravity and topography). 
The other models we are aware of, typically obtained either with a different methodology or 
(sometimes and) in different conditions (like spin coating e.g.), are different, see for instance [3H 
[511 H51 [5DI [55] . On the other hand, recent works in the mathematical literature also studied 
reduced models for thin layers of viscoelastic flows. For instance [71 [5] derive reduced models for 
the Oldroyd-B (OB) system of equations, where a purely viscous component div {'risD{u)) is added 
to the stress term in the right-hand side of (|2.2p in comparison with the UCM equations. But 
our project is different in essence from the thin layer models obtained for those viscoelastic flows 
without free surface and essentially driven by viscosity instead of gravity. Recall that here we 
focus on gravity-driven shallow regimes, and that is why we consider the UCM model in particular 
rather than the OB model (the viscosity only plays a minor role here). 

In Section [5] below, we recall the UCM model for viscoelastic fluids and some of its properties 
in the mathematical setting that is adequate to our model reduction. Then our new reduced model 
is derived in Section [3] under a given set of clear mathematical hypotheses (which we classically 
cannot embed into an existence theory for solutions to the non- reduced UCM system of equations). 
Section 2] is devoted to the study of some mathematical properties of our new reduced model. In 
Section [51 we provide numerical simulations in benchmark situations where shallow viscoelastic 
flows could be advantageously modelled by our new system of equations. Last, in Section \6\ a 
physical interpretation of situations modelled by our system of equations is given in conclusion, 
along with threads for next studies. 

2. Mathematical setting with the Upper-Convected Maxwell model for 

viscoelastic fluids 

The evolution for times t E [0, +oo) of the flow of a given portion of some viscoelastic fluid 
conflned in a moving domain Pt C M"* (d = 2 or 3) with piecewise smooth boundary OVt is governed 
by the following set of equations, the so-called Upper-Convected Maxwell (UCM) model [H[^l42j: 

(2.1) divM = in 2?t, 

(2.2) dtU+{u-V)u^-Vp + divT + f in Pj, 

(2.3) dtT + {u-V)T = {Vu)T + T{Vuf + \{r]pD{u)-T) in A, 

A 

where: 

• u : {t,x) e [0, +oo) X I>t i-> u{t, x) e W- is the velocity of the fluid, 

• D{u) : {t,x) e [0,+oo) X K> D{u){t,x) G R^'"', where M^'"' denotes symmetric real 
dxd matrices, is the rate-of-strain tensor linked to the fluid velocity u through the relation 

(2.4) £)(m) = i(VM + Vm^). 

• p : {t,x) ^ (0, +oo) X I?t p{t, a;) e M is the pressure. 



3 



• r : (t, a;) G [0, +00) x 2?t 1— > T{t, x) G R^^'' is the symmetric extra-stress tensor, 

• r]p, A > are physical parameters, respectively a viscosity only due to the presence of 
elastically deformable particles in the fluid, and a relaxation time corresponding to the 
intrinsic dynamics of the deformable particles, 

• f : {t,x) e [0, +00) xVt^ f{t, x) e R'' is a body force. 

Notice that we have assumed the fluid homogeneous (with constant mass density, hence normalized 
to one). We also refer to the Section [6] where more details about the UCM model are given along 
with a physical interpretation of our results. From now on, we assume translation symmetry 
(c? = 2), we endow with a cartesian frame (e^,, e^) such that / = —ge^ corresponds to gravity 
and we assume that Vt has the following geometry (in particular, surface folding like in the case 
of breaking waves is not possible) : 

(2.5) Vie [0, +00), X = ix,z) eVt <^ X e {0,L), < z -b{x) < h{t,x), 

where b{x) is the topography elevation and b{x) + h{t, x) is the free surface elevation of our thin 
layer of fluid. Note that the width h{t, x) is an unknown of the problem (it is a free boundary 
problem). We shall denote as ax (respectively a^) the component in direction (resp. e^) of any 
vector (that is a rank-1 tensor) variable a, and similarly the components of higher-rank tensors : 
(ixxyCixzi ■ ■ ■ We denote by n : x £ (0,L) — > n{x) the unit vector of the direction normal to the 
bottom and inward the fluid: 

(2.6) n 



\/l + {dxbr VlTTW 
We supply the UCM model with boundary conditions for all t E (0, +00): pure slip at bottom, 

(2.7) un = Q, forz = 6(a;), xe{0,L), 

(2.8) rn = {{rn) ■ n)n, for z = 6(a;), xe(0,i), 

kinematic condition at the free surface Nt + N ■ u ~ Q where {Nt , N) is the time-space normal, i.e. 

(2.9) dth + Uxdxib + h) ^ Uz, for z = b{x) + h{t, x), xe{0,L), 
no tension at the free surface, 

(2.10) {pi -T)-{-dxib + h), 1)^0, for z ^b{x) + h{t,x), xe{0,L), 

plus (for example) inflow/outflow boundary conditions or periodicity in x. We insist on (j2.8p 
without friction. Adding a friction term in (j2.8p would not yield the same result. Finally, the 
Cauchy problem is supplied with initial conditions 

(2.11) m(0, x) = u°{x), t(0, x) ^ T°{x), /i(0, x) = h°{x), 

assumed sufficiently smooth for a solution to exist. Note indeed that the existence theory for 
solutions to the UCM system ((^IM^^MO)) is still very hmited (see e.g. [551 IM]), like for non- 
Newtonian flows with a free surface (see e.g. [5^ [55] for the so-called Oldroyd-B model with a 
viscous term in (|2.2p ). 
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Last, we recall some essential features of the UCM model (j2.1H2.11"|) . Let cr : {t, x) E [0, +00) x 
Vt H> (T{t,x) € Rg^"^ be the symmetric conformation tensor linked to the symmetric extra-stress 
tensor t through the relation 

9\ 

(2.12) rr = I+—T, 

where I denotes the c?-dimensional identity tensor. The UCM model can be written using the 

Ih. , 

2\ 



variable cr instead of t. Indeed, diver replaces divr in (j2.2p . and (|2.3p should be replaced with 



(2.13) dt(T + {u-V)(T = {Vu)<T + <T{Vuf + \{I-<t) in 2?t . 

In addition, the following properties are easily derived following the same steps as in |16) for the 
Oldroyd-B model (except for the absence of the dissipative viscous term r]s\D{u)\'^). First, for 
physical reasons, cr should take only positive definite values (this is easily deduced when cr is 
interpreted as the Grammian matrix of stochastic processes, see [M] e.g. and Section W^ . The 
initial condition (j2.1ip should thus be chosen so that (T{t = 0) is positive definite. Provided 
the system (j2.1H2.11"|) has sufficiently smooth initial conditions and the velocity field u remains 
sufficiently smooth, cr indeed remains positive definite (see [TS] e.g.; where the viscosity rjs plays 
no role in the proof). Second, the system (|2.1H2.11"|) is endowed with an energy (the physical free 
energy) 

(2.14) F{u,t) = J^ (^^\u\^ + ^tr{cT-lna-I)-f-x^dx 
which, following Reynolds transport formula and |16) . is easily shown to decay as 

(2.15) |f(m, t) = ^ tr(<T + a-^ + 2I)dx . 

3. Formal derivation of a thin layer approximation 

Our goal is to derive a reduced model approximating (j2.1H2.11"|) in the thin layer regime /i <C L 
where L is a characteristic length of the flow. We follow the formal approach of [111 [T5] . Our main 
assumption is thus h/L = 0(e), thereby introducing an adimensional parameter e — ^ 0. In the 
following, we simply write 

(HI) /i = 0(e), 

where one should infer the unit (L) and the limit (e — > 0). This corresponds to a rescaling of 
the space coordinates with an aspect ratio e between the vertical and horizontal dimensions. We 
shall also implicitly use a characteristic time T over which the physical variables vary significantly. 
Then, our task can be formulated as: find a set of non-negative integers 

such that a closed system of equations for variables {u,p,f) approximating {u,p,t) holds and 
(3.1) {u-u,p~p,T-T) = 0{e^) 
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is a uniform approximation on Vt (where the powers are apphed componentwise). Note that 
in p.ip and all along the paper, O has to be understood componentwise in the unit corresponding 
to the variable. In particular, the previous assumption p.ip more explicitly signifies 



|(«-«),|J(p-p),|J(r-T)) =0(6^). 

We proceed heuristically, increasing little by little the degree of our assumptions on /. Hopefully, 
the reduced model found that way corresponds to a physically meaningful regime. 

We recall that another viewpoint is to find a closed system of equations for depth-averages of the 
main variables of the system (j2.1H2.11"|) . The link with our approach is as follows. The conservation 
of mass for an incompressible, inviscid fluid governed by (|2.ip within a control volume governed by 
the evolution of the free-surface height as given by the kinematic boundary condition (|2.9p and the 
boundary condition (|2.7p at the bottom reads as an evolution equation for the free-surface height 
h (using the Leibniz rule) where the depth-averaged velocity profile Ux '■— \ /d'^'' Uxdz enters, 
(3.2) 

pb+h / I'b+h \ 

'it, X e [0, +oo) X (0, i) = / {dxUx + dzUz)dz = dth + dx \ Uxdz 1 = dth + dx{hux) ■ 

The challenge in the derivation of a reduced model for thin layers is then to find a closure for the 
evolution of iix in terms of the variables {u,p, t). In particular, depth-averaging the equation for 
Ux with the boundary conditions (|2.9H2.8H2.10p gives, using again the Leibniz rule, 

/ I'b+h \ / .b+h \ 

(3.3) ^My Uxdzj+dxij (ul+p-Txx) dzj =[{Txx-p)dxb-Txz]\b, 

showing that a typical problem is to write an approximation for jj^^'^ u^. and for the source term 
in the right-hand-side of (j3.3p in terms of {u,p, t). 

We now give the detailed system of equations (|2.1H2.11"|) in the 2-d geometry of interest: 



(3.4a) 
(3.4b) 
(3.4c) 


dxUx + dzUz = 0, 

dtUx + UxdxUx + UzdzUx = 

dtUz + UxdxUz + UzdzUz = 


-dxP + dxTxx 
-d,p + dxTxz 


+ dzTxz, 

+ dzTzz - g, 




(3.4d) 


^t^xx ~l~ ^x^x^xx H" '^z^z'^xx 


— i^'^^x'^x^'^xx 


+ {2dzUx)Txz 


i ^ ^X ^X ^ ' XX 1 


(3.4e) 


dtTzz + UxdxTzz + uATzz 


= {2dxUz)Txz - 


f {2dzUz)Tzz ■ 


11 1 

4- -^dzUz - jTzz, 


(3.4f) 


dfTxz + UxdxTxz + UzdzTxz 


= {dxUz)Txx 4 


- {dzUx)Tzz + 


^(dzUx + dxUz) - -^Txz, 



6 



where we have used p.4a|) to simphfy (|3.4fp . The boundary conditions (|2.7p . (|2.8p and (|2.10p write: 

(3.5a) Uz = {dxh)ux at z = 6 , 

(3.5b) - idxb)Ta:x + T^z = -dxb{-{dxh)Txz + Tz^ at z = 6 , 

(3.5c) - dx{b + h){p - Tx.j,) - T^z ^ Q &tz^b + h, 

(3.5d) dx{b + h)Txz + {p - Tzz) ^ a.t z = b + h, 

while the kinematic condition (|2.9p . foUowing p.2p . writes 

(3.6) dth + d, ^ ^ w,dz^ ==0. 

We first simphfy the derivation of a thin layer regime by assuming that the tangent of the angle 
between n and is uniformly small 

(H2) d^b ^ 0(e) as e ^ , 

hence only smooth topographies with small slopes are treated here. This restriction could prob- 
ably be alleviated following the ideas exposed in [TS], though at the price of complications that 
seem unnecessary for a first presentation of our reduced model. On the contrary, the following 
assumptions are essential: 

(H3) ?7p - 0(e), A -0(1). 
(As explained previously, we recall that the assumptions (H3) hold in the unit of the variable, 
which is here L^/T and T respectively.) As usual in Saint Venant models for avalanche flows, we 
are looking for solutions without small scale in t and x (thus only with scales T and L), but with 
scale of order e in z (in fact, ei), which can be written formally as 

(3.7) 9^^=0(1), 9. -0(1), 9, = 0(l/e) 

in the respective units l/T, 1/L, 1/L. From now on, for the sake of simplicity, we shall not write 
explicitly the units as functions of T and L wherever they come into play. 

We are looking for solutions with bounded velocity u with bounded gradient Vtt. Thus accord- 
ing to (|3.7p and to p.5ap . we are led to the following assumptions on the orders of magnitude 

(H4) = 0(1), Uz = 0(e), dzu, = 0(1), as e ^ 0. 

(A typical profile for reads A{t/T, x/L) + zB{t/T, x/L), with any dimensional functions A and 
B of the adimensional variables t/T and x/L.) 

According to (j2.3p . a typical value for t is ripD{u). Thus we assume accordingly that 

(H5) T ^ 0(e) as e ^ 0. 

We deduce from above that there exists some function u"(i, a;) depending only on {t,x) such that 

(3.8) Ua,{t,x,z)=ul{t,x)+0{€). 

Then, following the classical procedure |26l[Tn[T51l37) . we find the following successive implications. 



7 



i) From the equation (|3.4cp on the vertical velocity Uz, we get by neglecting terms in 0(e) 

(3.9) dzp = dzTzz -g + 0{e) . 

Hence d^p — 0(1), and the boundary condition p.5d[) gives that p = 0(e), indeed 

(3.10) p^Tzz+9{b + h^ z) + 0{e^). 

ii) Next, from the equation (|3.4b[) on the horizontal velocity we get 

(3.11) dtul + uld.,ul^dzT.,z + 0{e). 

The boundary condition (j3.5bp gives Txz\z=b ~ O(e^), thus with p. lip it yields 

(3.12) T^z = [dtul + uld^uDiz -b) + 0{^) . 

In addition the boundary condition p.5cp implies that Txz\z=b+h ~ 0{e^). We conclude 
therefore that 

(3.13) dtul + uld^ul = 0(e), T^z = 0(e2). 

iii) The previous result combined with the equation p.4f|) on t^z implies dzUx = 0(e), hence 

(3.14) Ux{t,x,z) = ul{t,x) + 0{e^). 



This "motion by slices" property is stronger than the original one 
iv) Using ((XTi)) and (|XTn)) in (|T4b| improves ([XTT|) to 

(3.15) dtul + uldxul = dx{Txx - Tzz - gib + h)) + dzTxz + O(e') , 

which gives, with the boundary condition (j3.5bp [t^z ~ dxb{Txx — Tzz)] \z=b — 0(e'^ 

Txz = [dxb{Txx - Tzz)] \z=b " / dx{Txx ~ Tzz) dz 



(3.16) 

+ {dtul + uldxul + gdx{b + h)) (z - 6) + 0(e3). 

But according to p.5cp combined with p.lOp . one has [txz — dx{b + h){Txx ~ Tzz)] |z=6+/i = 
0(e3), thus with ((XT5|) 



(3.17) 



\dx{b ^ h){Txx - Tzz)\\z=b^h - \ dx(Txx - Tzz) dz 



b+h 



+ {dtul + uldxul + gdx{b + h)){z^b~h) + Oie^). 
Therefore, the difference of p.l6p and p.l7p yields 

(3.18) {dtul + uldxul + gdxib + h))h^ dx (^J^ (txx - Tzz) dz^ + 0(e3) . 

We note that Txz is then given by p.l6p or (j3.17p as a function of m° and {txx — Tzz), and 
the evolution equation p.lSp for u° is exactly the one that one would have obtained after 
integrating p.lSp in the direction and using the boundary conditions p.Sbp and p.Scp 
combined with p.lOp . It can also be obtained from 
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v) The result (j3.14p with the incompressibihty condition (j3.4ap and the impermeability condi- 
tion p.5ap at the bottom also allows to compute the vertical component of the velocity 



(3.19) 



which is of course consistent with our hypotheses about u, = 0(e). 
vi) Collecting all the previous results, (|3.4dp and p.4ep up to O(e^) give 
(3.20) 



dtTxx + u^dxTxx + {{dxh)ul - {z- b)dxu°) d^r^x = 2{dxul)T^x + 



A 



dtT. 



u°5,r,, + {{dxb)ul -{z- b)dxul) a.r,, = -2(9.^°)r,, - ^^^-"^ + + Qie'), 



which closes the system of equations for the reduced model, 
vii) The previous results which give Txz at order 0{e^), that is p.l6p or p.l7p . are consistent with 
the equation (j3.4fp for t^z at order O(e^), from which one could next obtain an approximation 
for dzUx up to O(e^), that is 



(3.21) dtT.z + uldxTxz + {{dxb)ul + (Z - b)dxu'l^)dzTxz + \t,z 



Z \ I zz 



Ik] 

2XJ 



= dx{idxb)ul{z-b) + dxul) (txx . 

with and Tzz given up to order O(e^) by p.20p . This procedure fixes the next term in 
the expansion p.l4p . Note in particular that we do not have Ux{t,x,z) = u'^{t,x) + 0{e^) 
(dependence on the vertical coordinate subsists at order e^). 

To sum up, dropping e, we have obtained a closed system of equations 

( dth + d,ihul)^0, 



(3.22) 



at(K) + a. U(<)' + .9Y + 



b+h 



dtTxx + uldxTxx + {{dxb)ul ~ {z~ b)dxul) d, 

dtTzz + uldxTzz + {{dxb)ul - (Z - b)dxul) dz 



-g{dxb)h, 



Txx = 2{dxul)T,, + '^dxu'i ~ 
Tzz = -2idxul)Tzz - '-^dxul 



A 



which allows to compute consistently uniform asymptotic approximations of (m^, Uz,p, Txx-, Tzz, Txz) 
as variables of order 0(e^°'^'^'^'^'^-'), up to errors in 0(£(2:3,2,2,2,3)-j^ These correspond to approxi- 
mations of dUiD-dini) up to 0(e(2.2,l,2,2,3,3,3,3,2,3))_ 

In (j3.22p . b depends only on x, h and it*J depend on (i,a;), while t^-x and Tzz depend on 



(t, X, z). However, observe that the momentum conservation equation invokes only ^^^^ T^xdz and 



Ib^^ Tz^dz, which do not depend on z. Now, using Leibniz rule and the boundary conditions, 



it is possible to get equations for J^^ ^ r^xdz and Tzzdz (integrating those for t^x and r^^. 
see (|4.4p below in Section U) and form a closed system with the equations for the momentum and 
mass conservation. Another equivalent way to derive the same closed system of equations is to 
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assume that Txx and Tzz are independent of z (at least at first-order in e). In the rest of this 
paper, we shall mainly be concerned with that simplified system of equations, whose mathematical 
properties are easier to study. 

4. The new reduced model and its mathematical properties 

The reduced model p.22p is endowed with an energy equation similar to the one for the full 
UCM model. Obviously, the whole system of equations for r in the reduced model rewrite with 



2A. 

we consider only the diagonal part 

(4.1) ^,^hy.. = i + -r.. 

The two last equations of (|3.22p yield 



the entries of the conformation tensor a = / + — r. However, since it is diagonal at leading order, 



1 I 2A 

O zz ' 1 n~ 'zz , 
'IP / 



(4.2) 



dt<7xx + u°dx<Txx + {{dxb)ul - (z - b)dxu°) d^axx = 2{dxul)axx - ^{cTxx - 1), 
dta.. + uldxu,, + {{dxh)ul - (z - b)dxul) d,a,, = -2{dxul)a,, ~ ^{a,, ~ 1). 



These equations imply that cixx and cizz remain positive if they are initially. Then, we compute 
(4.3) 

dt + uldx + {{dxb)ul - (Z - b)dxul) dz) {]-Txx - 7^ ln(l + —Txx)) = {dxul)Txx - 

[dt + uldx + {idxb)ul -{z- b)dxul) dz) {\tzz - ^ ln(l + ^r,,)) = -(9.«°)r,, " ^ J- 

In order to compute the integral of (|4.3p with respect to 2, we notice the following formula for 
any function ipit, x, z) (a combination of the Leibniz rule with boundary conditions at z = b and 

z = b + h), 

rb+h 

[dt + uldx + {{dxb)ul - (z - b)dxul) dz) ^ dz 



b. 



b+h 



b 



(dt^ + dxiul^)+dz[{idxb)ul -{Z- b)dxul) dz 
Hi pb+h 



(^■"^^ = dt I ip dz - ipb+hdth + dx ulip dz - {ul(p)b+hdxib + h) + {ulip)bdxb 



'(^{dxb)ul - hdxul)tpb+h - idxb)ul(pb 



fb+h / fb+h 

dt I tpdz + dx ( u° J ipdz 



Therefore, summing up the two equations of (j4.3p and integrating in z gives 
(4.5) 



b+h -, i-b+h / 2 2 

{dxuD / (r,, -Tzz)dz~- {^ + -^]dz. 

Jb Vp Jb 
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Moreover, the classical computation of energy for the Saint Venant model gives 



(4.6) dt(^i 



-g-+gbh) + d,[{h^-^ + gh' + gbh)ul 



b+h 



{Tzz ~ T^.x)dz ^ 0. 



Adding up and yields 



(4J) 






+ 9 


+dx (J^h 


{<? 




2 




rb+h 
/ i- 



4A2 



gh^ + gbh + ^ / tr(ffO - In ^7° - J) rfz + ^ / 
4A Ji, 2a Ji, 

rb+h 

/ tr(cr" + [crO]-i -27)rfz. 

Jb 



Therefore, we get an exact energy identity for solutions to the reduced model p.22p . Note that to 
discriminate between possibly many discontinuous solutions (generalized solutions in a sense to be 
defined, see below the discussion on the conservative formulation), we would naturally require an 
inequality in (j4.7p instead of an equality. 

In the case of t^x and Tzz independent of z, everything becomes more explicit. Using the 
variables Gxx — i ' 
reduced model then writes 



1 + —Txx and Uzz = 1 + Tzz (also clearly independent of z), the simplified 



(4.8) 



' dth + dx{hu°) =0, 
dtihul) + dx {h{ulf + + ^/^(fT.. - ^xx 

1 - (^xx 

X ' 

1 — (T,, 



-ghd^b, 



dto-xx + u°j,dx(7xx - 2axxdxul = 
dt<Jzz + u^xdxCTzz + 2azzdxul = 



A 



while the energy inequality becomes (cr'^ is defined in (j4.1l) ) 



dt h 



+ gbh + ^htria" - ln<TO - I 
2 ^ 2 4A ^ 



(4.9) 



,0\2 



2 , gh^ + gbh + '^h trier" - In^r" - /) + ^h{azz ^a^x)] w° 



<-^htr{(T" + [<T"]-'~2I). 
- 4A2 

In (|4.8p and (|4.9p . 6 is a function of x and h, axx, (^zz depend on (i, x), with /i > 0, cr^x > 0, 
Czz > 0. From now on, we shall only deal with the simplified reduced model (j4.8p . 

The inequality (j4.9p (instead of equality) for possibly discontinuous solutions rules out general- 
ized solutions for which the dissipation - already present in our model ! - is physically not enough 
(see also |16) where a similar numerical "entropy" condition is used to build stable finite-element 
schemes for the viscous UCM model, namely the so-called Oldroyd-B model). 
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Remark 1 (Limit cases). For the system {4-^ ! two interesting regimes are important to mention. 
The first is the standard Saint Venant regime, for which one takes rjp/\ = 0. It is obtained in the 
limit r]p ^ for fixed X (as opposed to the limit A — (X) for fixed rjp, some kind of "High-Weissenherg 
limit" |42| which is problematic, as we will see in the numerical experiments) . The second regime is 
obtained in the "Low-Weissenberg limit" X — > 0, for fixed rjp. Assuming {l~(7xx)/X and (l'~azz)/X 
remain bounded, the system rewritten with t^x o,nd t^z gives the viscous Saint Venant system 



( dth^-d^{hul) = {), 



(4.10) 



dt{hul) + d., ( h{ulf+g^ ~ 2r]phd^ul ) - -ghd^h, 



with the energy inequality 
(4.11) 

Remark 2 (Steady states). The source terms (1 ^ a^^) / X and (1 ^ azz)/ X in ^.8^ are responsible 
for the right-hand side that dissipates energy in ^■9\ l. This dissipation has the consequence that 
steady states are possible only if 

(4.12) tr(cr" + [cr"]"^ - 21) = 0, i.e. x = , 

which implies that steady solutions to ^.8^ identify with the steady solutions at rest to the standard 
Saint Venant model: u'^ = 0, h + b = est, a^x — f^zz ~ 1- 

Remark 3 (Conservativity). The reduced model {J^.S^ is a first-order quasilinear system with 
source, but not written in conservative form because of the stress equations on a^x md Ozz- Indeed, 
one can put them in conservative form as follows. 



(4.13) 



dt {{axx)-'/') + {{-xx)-'/'ul) = -a-^/^ 
d,{{azzf'')+d.^{{azzY''^:) = 



1 — CFxx 



2X 



-111 1 ~ 



a 

2A 



However, these conservative equations do not help since they are physically irrelevant. Moreover, 

1/2 

the physical energy of {4-^^ *s ^oi convex with respect to these conservative variables Oxx and 

1 /2 

azz ■ As a matter of fact, one can show that the energy, that is 

(4.14) E = + ffy + + (^^- + - H'^^x'Jzz) - 2) , 
cannot be convex with respect to any set of conservative variables of the form 

(4.15) (^h, hul hw-' 1^-^^ , h,-' 1^-^^ ^ , 

where zu,i; are smooth functions standing for general changes of variables, see AvvendixTM 
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Nevertheless, the system (j4.8p can be written in the following canonical form, strongly reminis- 
cent of the gas dynamics system, 

( dth + d^ihu",) = 0, 

(4.16) { ^'('^"S) + ('^("S)' + PiK s)) = -ghd^b, 

1 



with 
(4.17) 



dts + uldxs = ^S{h,s) 



' -1/2 1/2^ 



(4.18) 



-3/2 -1/2 
'^XX X ^zz 

-(1 - Oxx), 



2h 



2h 



(1 



(4.19) 

One can compute 
(4.20) 



dP 

dh 



— ) =9h 



P{h,s) 



2A 



Ha^z 



2a, 



from which we conclude that for smooth b, the system (j4.16p is hyperbolic with eigenvalues 
(4.21) Xi^ul- 



gh 



>(3a,, + a,,), A2 = ul A3 - w° + Jgh + ^(3a, 



+ CTxx 



the second having double multiplicity. One can check that A2 is linearly degenerate, while Ai and 
A3 are genuinely nonlinear (this follows from computations similar to [27j Example 2.4 p. 45] and 
the first line of 

From the particular formulation (|4.16p . one sees that the jump conditions for a 2— contact 
discontinuity are that and P do not jump (as weak 2-Riemann invariants). However, jump 
conditions across 1— and 3— shocks need to be chosen in order to determine weak discontinuous 
solutions in a unique way. 

A possible choice of jump conditions is, as explained in Remark [31 to take the conservative for- 
mulation (|4.13p (or equivalently a conservative formulation related to the variables (|4.15p , leading 
to the condition that s does not jump through 1— and 3— shocks). This formulation gives unphys- 
ical conservations and nonconvex energy (which could produce numerical under/overshoots), and 
we shall not make this choice. 

Our choice of jump conditions will be rather imposed indirectly by numerical considerations, 
via the choice of a set of pseudo-conservative variables, i.e. variables for which we shall write 
discrete flux difference equations. Solving nonconservative systems leads in general to convergence 
to unexpected solutions, as explained in |18j . With a pragmatical point of view, we nevertheless 
choose the pseudo-conservative variables as 



(4.22) 



q= (gi, 92, 93, 94)^ := {h,hul,haxx,h(T,z)' 
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In other words, we consider the formal system 

dt{hul) + a, {h{ulf +9Y + ^^(^- - '^--)) = '9hd, 

h — hoxx 



(4.23) 



dtiha^x) + dxiha^xul) - 2haxxdxU° = 
dt{hazz) + dx{hazzU°) + 2hazzdxU° = 



A 

h - hffzz 
A ■ 



The choice of these pseudo-conservative variables is good for at least two reasons: 

• these variables are physically relevant, 

• the energy E in (|4.14p is convex with respect to them (see Appendix |X| . 

The second point will make it easier to build a discrete scheme that is energy satisfying (in the 
sense of the energy inequality (|4.9p ). while preserving the convex (in the variable q) set 

(4.24) U = {h> 0, <Jxx > 0, <Jzz > 0} , 

which is here the physical invariant domain where the energy inequality (|4.9p makes sense. Note 
that our system is of the form considered in (see also Remark U]). 

Let us mention that for the viscous UCM model, namely the Oldroyd-B model, various numer- 
ical techniques are proposed in [IHl 13 for the preservation of the positive-definiteness of a 
non-necessarily diagonal tensor cr in the context of finite-element discretizations. 

5. Finite volume method and numerical results 



In this section we describe a finite volume approximation of (j4.23p . The approximation of the 
full system is achieved by a fractional step approach, discretizing successively the system (|4.23p 
without the relaxation source terms in 1/A on the right-hand side of the two stress equations, 
and these relaxation terms alone. The topographic source term hdxb is treated by the hydrostatic 
reconstruction method of |3] in Subsection 15.41 This approach ensures that the whole scheme is 
well-balanced with respect to the steady states of Remark [21 because the relaxation terms vanish 
for these solutions. 

The integration of relaxation source terms is performed by a time-implicit cell-centered formula. 
Note that then the scheme is not asymptotic preserving with respect to the viscous Saint Venant 
asymptotic regime A — >■ of Remark [1] for this one would need a more complex treatment of these 
relaxation terms. 

Let us now concentrate on the resolution of the system (j4.23p without any source, i.e. 

' dth + dx{hul) = 0, 
dtihu^'J + dx {Hulf + P) = 0, 
dt{haxx) + dx{haxxul) - 2haxxdxui = 0, 
^ dtihozz) + dx{hazzu°) + 2h(Tzzdxul = 0, 



(5.1) 
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with 



(5,2) 




and the energy inequahty 



(5.3) 




ZZ 



, -ln(a,,a,3)-2)+p)<j < 0. 



A finite volume scheme for the quasihnear system (|5.ip - (|5.2p can be classically built following 
Godunov's approach, considering piecewise constant approximations of q = {h^hvP^,haxx,hazz), 
and invoking an approximate Riemann solver at the interface between two cells. 

5.1. Approximate Riemann solver. In order to get an approximate Riemann solver for (|5.ip . 
we use the standard relaxation approach, as described in |12| . It naturally handles the energy 
inequality (j5.3p , and also preserves the invariant domain (j4.24p . 

Because of the canonical form of (|5.ip . which is (j4.16p without source, i.e. 



we have a formal analogy with the system of full gas dynamics equations. Therefore, we follow 
the usual Suliciu relaxation approach that is described in [T^]. We introduce a new variable tt, 
the relaxed pressure, and a variable c > intended to parametrize the speeds. Then we solve the 
system 



(5.4) 




with 



(5.5) 




(5.6) 



< 



' dth + dx{hul)={), 
dt{hul)+dx{h{ulf +Tr) = Q, 
dtihn/c') + dxihnul/c' + ««) = 0, 



dtc + u°dxC = 0, 
dfS + u'^dxS = 0. 
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This quasilinear system has the property of having a quasi diagonal form 




dtc + u^.dxC = , 
dts + u^dxS — . 



One deduces its eigenvalues, which are vP^ — c//i. w^J + c//i, and u*J with multiplicity 4. One checks 
easily that the system is hyperbolic, with all eigenvalues linearly degenerate. As a consequence, 
Rankine-Hugoniot conditions are well-defined (the weak Riemann invariants do not jump through 
the associated discontinuity), and are equivalent to any conservative formulation. We notice that 
with the relation (jS.Sp the equation on s in (|5.6p can be transformed back to 



The approximate Riemann solver can be defined as follows, starting from left and right values of 
h, hu^, haxxi huzz at an interface : 

• Solve the Riemann problem for (|5.6p with initial data completed by the relations 



and with suitable values of c; and that will be discussed below. 
• Retain in the solution only the variables h, hu^, ha^x, hcTzz- The result is a vector called 
R{x/t,qi,qr). 

Note that this approximate Riemann solver R(x/t, qi, qr) has the property to give the exact solution 
for an isolated contact discontinuity (i.e. when the initial data is such that and P are constant), 
because in this case the solution to (|5.6p is the solution to (jS.ip completed with tt = P{h, s). 

Then, the numerical scheme is defined as follows. We consider a mesh of cells {Xi_ii2,Xi^i/2), 
I £ Z, of length Axi = Xi+1/2 — 2:1-1/2) discrete times t„ with — tn = At, and cell values 

approximating the average of q over the cell i at time tn- We can then define an approximate 
solution q'^PP''{t,x) for tn < t < tn+i and .t £ R by 



(5.8) 



dtihcFxx) + dx{haxxUx) - 2h(JxxdxU°^ = 0, 
dt{hazz) + dxiha^zUx) + 2hazzdxU% = 0. 



(5.9) 



TT; = P{hl, {(Txx)l, {(yzz)l), 



TTr = P{hr, {cTxx),-, {(Jzz)r), 



(5.10) 




where Xi 



(2:^-1/2 + 2;i-|-i/2)/2. This definition is coherent under a half CFL condition, formulated 



as 



(5.11) 
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The new values at time are finally defined by 

(5.12) gf+i = r^^'^' q^PP^itr.+, ~ 0, x) dx. 

Notice that this is only in this averaging procedure that the choice of the pseudo-conservative 
variable q is involved. We can follow the computations of Section 2.3 in |T2], the only difference 
being that here the system is nonconservative. We deduce that 

(5.13) g,r+' = - ^ imi^cQ+i) M<l7-i.<l7)) , 
where 

^MuQr) ^ F{qi) - I (R{^,qi,qr) - qi)d^, 

(5.14) J-^^ ^ 

J^riqi,qr) = F{qr) + J i^R{^ , qi , qr) - qrj d£„ 

and the pseudo-conservative fiux is 

(5.15) F{q) = {hul, hiulf + P, ha,.,ul, ha,,ul). 

In (|5.15p . the two last components are chosen arbitrarily, since anyway the contributions of F in 
(|5.13p cancel out. 

Since the two first components of the system (j5.6p are conservative, the classical computations 
in this context give that for these two components, the left and right numerical fluxes of (|5.14|) are 
equal and indeed take the value of the flux of (|5.6p . i.e. hvP^ and hiu^Y + tt, at x/t = 0. 

We can notice that while solving the relaxation system (j5.6p . the variables /i, Sxx and Szz remain 
positive if they are initially (indeed this is subordinate to the existence of a solution with positive 
/i, which is seen below via explicit formulas and under suitable choice for ci, c^). By the relation 
(|5.5p this is also the case for and Uzz- Therefore, the invariant domain lA in (|4.24p is preserved 
by the numerical scheme (j5.13p . this follows from the average formula (j5.12p and the fact that U 
is convex (in the variable q). 

Remark 4. The above seheme satisfies the maximum principle on the variable Sxx, o-nd the mini- 
mum principle on the variable Szz- This means that if initially one has Sxx < k for some constant 
fc > (respectively Szz ^k), then it remains true for all times. 

This can be seen by observing that the set where Sxx < k (respectively Szz > k) is convex in 
the variable q, because according to (j5.5p . (j4.22p . it can be written as qiq^ > (respectively 
k^qf — qi < 0)- Then, s is just transported during the resolution of (j5.6p . while the averaging 
procedure (j5.12p preserves the convex sets. Another proof is to write a discrete entropy inequality 
for an entropy h(j)(sxx), which is convex «/ < < Sxx4'" ; take for example 4>{sxx) — max(0, Sxx ~ 
A;)^/2 (respectively for an entropy h(j)(szz), which is convex if < ^<f)' < iszz4>" , take for example 
(j)(szz) = k^^^'^Szz — f'Szz'^ + ^k'^^^ for Szz < k, (j){szz) ~ for Szz > k). We shall not write down 
the details of this alternative proof. 
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5.2. Energy inequality. We define in a similar way the left and right numerical energy fluxes 

gi{qi^q,)^G{qi)- f ( E {R{t qi , Qr)) - E (qi)) , 

(5.16) J-JS.) \ 
Gr{quqr)^G{qr) + i^E{R{tqi,qr))-E{qr)jdC, 

where E is the energy of (|4.14p without the topographic term gbh, and 

(5.17) G = (i? + P)«° 

is the energy flux. We have from [T2] that a sufficient condition for the scheme to be energy 
satisfying is that 

(5.18) gr{qi,qr)-Gl{qi,qr) <0- 

When this is satisfied, because of the convexity of E with respect to q one has the discrete energy 
inequality 

(5.19) E{q-+') - EicS) + ^(c;(g»,g7^^) - g{q-^,,q-)) < 0, 

where the numerical energy fiux g{qi, qr) is any function satisfying Griqi, qr) < Q{qi, qr) < Gi{qi, qr)- 
In order to analyze the condition (|5.18p . let us introduce the internal energy e{q) > by 

h 77 

(5.20) ^ ^ ^ 2 4A ^'^'"'^ ^ '^^^ ^ ^T^'^i'^xxCTzz) - 2) , 
so that 

(5.21) E ^h{ulf 12 + he, 

and {dhe)\a = P/h?. Then, while solving the relaxation system (|5.6p . we solve simultaneously the 
equation for a new variable e. 

(5.22) dt{e- 7rV2c2) + uld^{e^ tt^ /2c^) ^ 0, 

where e has left and right initial data e{qi) and e{qr). The reason for writing (j5.22p is that 
combining it with (j5.6p yields 



(5.23) dt(h{ulY/2 + h?^ +c),((/i(u°)V2 + /ie + 7r)uO) =0. 
Define now 

(5.24) g{quqr)=({h{ulfl2 + he + iTy\ . 

\ / x/t—0 

Lemma 1. If for all values of x/t the solution to (|5.6|) , (j5.22p satisfies 

(5.25) e>e(g), 

where here q = R{x/t,qi,qr), then Griqhqr) < Q{quqr) < Giiqi,qr) o,nd the discrete energy in- 
equality (j5.19p holds. 
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Proof. Since (|5.23p is a conservative equation, one has 

g{qi,qr)^G{qi)- f ({h{ulf/2 + he){0-E{qi))d^ 

(5.26) J-SS; 

= G{qr) + [{h{ulf/2 + he) {£) - E{qr))dC 

Therefore, comparing to (|5.16p . we see that in order to get the resuh it is enough that for all ^ 

(5.27) E{R{^,qi,qr)) < {h{u'i.f/2 + he) (0, 

which is dOSl). □ 

In order to go further, we fix the following notation: in the solution to the Riemann problem 
for (|5.6p . there are three waves and two intermediate states, denoted respectively by indices Z, * 
and r, *. Then we have the following sufficient subcharacteristic condition (recall that dhP\s is 
given by (g^Ol))- 

Lemma 2. If ci, Cr are chosen such that the heights h^ , h* are positive and satisfy 

yhe[hi,ht] h^dhPUh,si)<cf, 

(5-28) 

V/i e h^dhP\sih,Sr)<cl, 
then (j5.25p holds and thus the discrete energy inequality (j5.I9p is valid. 

Proof. The arguments of decomposition in elementary dissipation terms along the waves used in 
Lemma 2.20 in [T^] can be checked to apply without modification. □ 

Lemma 3. Denote 



(5.29) Pi = P{husi), Pr = P{hr, s,), ai = ^dhP\sihi, si), a, = ^/KpUK^, 
and define the relaxation speeds ci, Cr by 



= a; + 2 maxf 0, u° -° 



max 



(5.30) 

^ = 0^ + 2 max[o,u!^_, 



(O, Pr - Pi) 

hiai + hrar 
ma.x(Q, Pi - Pr 

hiQl + hrUr 



Then the positivity and subcharacteristic conditions of Lemma [H are satisfied, and the discrete 
energy inequality (|5.19p holds. 

Proof From (|4.20p and (|5.5p we have 
(5.31) dhP\s^gh+'^(3{hs ^2 



2AV ' [hs^xY 
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Denoting ip{h, s) = hy/dhP\s, we compute 

h f , Vp 



2VKPrs\ 2AV - h^sl. 



(5.32) = / , (2gh + ^ (3{hs,,f + jr^) +gh+^( &{h 



1 ft^u , a'lPfu. ^2 



3gh + 6^{hs,,y 



2^Kp\s ^ ' ^ 
Therefore, we deduce that ip satisfies 

dhf\s > 0, 

(5.33) 'fiih, s) — 7> oo as /i — > oo, 

Following [Proposition 3.2] [13] with a = 2, we get the result. □ 

Remark 5 (Bounds on the propagation speeds). Lemma \^ is also valid with the formulas of 
[Proposition 2.18] |12| instead of (j5.30p . Here we prefer (|5.30p because in the context of possibly 
negative pressure P these formulas ensure the following estimate on the propagation speeds: 

(5.34) max (^^, ^) < C + \ul,\ + ai + a,) , 

with C an absolute constant. This follows from the property that \P\ < hdhP\s, which is seen on 

5.3. Numerical fluxes and CFL condition. The Riemann problem for the relaxation system 
(|5.6p . (|5.22p has to be solved with initial data qi, qr completed with (j5.9p . the relation (j5.5p . 
ei = e{qi) = ei, e,. = e(g,.) = e^, and (j5.29p . (j5.30p . The explicit solution is given, according to 
|12| . by the following formulae, ft has three waves speeds Si < S2 < S3, 

(5.35) Si = , - q//i;, S2=U°^^, S3 = + Cr//lr, 

and the variables take the value "1" for x/t < Si, "1*" for Si < x/t < S2, "r*" for S2 < x/t < S3, 
"r" for S3 < x/t. The "1*" and "r*" values are given by 

_ _ '^'"2,; + Crw2,r + TT/ - TTr * _ * _ C,.7rj + QTT^ - ClCr{u°^^ - U° ,) 

Ci + Or Ci + Cr 

hi Cl{ci+Cr) ' h* hr Cr{ci+Cr) 

(5.37) Ci = Q, Cj. ~ Cr, Si ~ Si, Sj. = Sr, 

(5.38) cTxx.i = '^xx,i ( , a*^ r = f^xx.r (^j^^ , cr*^ , = cr^^,/ ^-^^ 



(5.36) 




("°)r 


= iO: 




1 







(5.39) =ei 



2cf ' 2cf ' 2c2 ' 2c2 
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(5.42) 



Then we need to compute the left /right numerical fluxes (|5.14p that are involved in the update 
formula (j5.13p . Since the h and hu'^ components in (|5.6p are conservative, classical computations 
give the associated numerical fluxes, and we have 

(5.40) -F, = J-''"", J';'"""), J'r = J''™",J^;"'-%J^;''""), 

where the conservative part involves the Riemann solution evaluated at x/t — 0, 
(5.4f ) T'^ = (K),/,^o, = (Kuir + 7r),/,^o- 

More explicitly, (|5.4f p yields that the quantities between parenthese are evaluated at "1" if Ei > 0, 
at "1*" if Si < < S2, at "r*" if S2 < < E3, and at "r" if S3 < 0. As usual there is no ambiguity 
in the resulting value when equality occurs in these conditions. The numerical energy flux (j5.24p 
involved in (j5.19p can be computed in the same way. 

We complete these formulas by computing the left /right numerical fluxes for the variables haxx , 
ha^z from (|5.14p . 

jr^'iT,, ^ {h(JxxU^x)i +i'rL\ii{Q,Y.i)(^{h(Jxx)*i - [ha^x)}^ 

+ miYi{0, 1:2) [{ha xx)*r - {haxx)fj + min(0, S3) ^(/itTa;^:)^ - {haxx)*^, 
(5.43) 

J-^^-- = {haxxu% - max(0,Si)((/ia,,),* - {haxx)i) 

- max(0, 'T.2)[{haxx)*r - {hiJxx)Vj ~ niax(0, S3) [{haxx)r ~ {hcjxx)*^ , 

the huzz fluxes being computed with the same formulas, replacing "xx" by "zz". 
The maximal propagation speed is then 

(5.44) ^((7,,(?.)=max(|Si|,|S2US3|), 
and the CFL condition (|5.1ip becomes 

(5.45) AtA{q,,qi+i) < i min(Aa:^, Aa;.,+i). 
Not that with ((OI)) and we get 

(5.46) A{qi,qr) < C {\ul i\ + \ul J +ai+ a^) 

with C an absolute constant, bounding the propagation speeds of the approximate Riemann solver 
whenever the left and right true speeds remain bounded. This property is more general than the 
possibility of treating data with vacuum considered in |12j . 
We have obtained finally the following theorem. 

Theorem 1. Consider the system (|5.ip with the pressure law (|5.2p . and denote the pseudo- 
conservative variable by q = {h,hu'^,haxx,hazz)- Under the CFL condition (|5.45p . the scheme 
(j5.13p with the numerical fluxes J^i{qi,qr), J^riqi,Qr) defined above via (|5.40p . and with the choice 
of the speeds (|5.29p . (j5.30p . satisfies the following properties, 
(i) It is consistent with (j5.ip - (|5.2p for smooth solutions, 
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(5.47) 



(ii) It keeps the positivity of h, a^x, ^^zz, 

(Hi) It is conservative in the variables h and hu'^, 

(iv) It satisfies the discrete energy inequality (|5.19p . 

(v) It satisfies the maximum principle on the variable Sxx, and the minimum principle on the 
variable Szz, 

(vi) Steady contact discontinuities where = 0, P = est are exactly resolved, 

(vii) Data with bounded propagation speeds give finite numerical propagation speed. 

(via) The numerical viscosity is sharp, in the sense that the propagation speeds of the approxi- 
mate Riemann solver tend to the exact propagation speeds when the left and right states qi, qr tend 
to a common value. 

5.4. Topography treatment. Consider now our system (|4.23p with topography, but without the 
relaxation source terms, i.e. 

( dth + dxihul) = 0, 

dt{hul) + dx (^hiuiy + ffy + ^^(^- - '^-)) - -ghdxb, 
dt{haxx) + dx{haxxul) - 2haxxdxul = 0, 
_ dt{hazz) + dxihazzu'i) + 2hazzdxul = 0. 

With respect to the previous sections, the term —ghdxb has been put back, where the topography 
is a given function b{x). For (|5.47[) . the energy inequahty (14. 9p is modified only by the fact that 
there is no right-hand side. Thus it can be written 

(5.48) dtE + dxG<0, 
with 

(5.49) E{q,b) ^ E{q) + ghb, G{q,b) = G{q) + ghbu^, 

where E and G are given by (|5.2ip . (|5.20p , (|5.17p . Recall that the steady states at rest of Remark 
[5] are defined by 

(5.50) M° = 0, h + b = cst, axx=(^zz = ^- 
Our scheme for (|5.47p is written as 

(5.51) q-+' = - ^ {FM7,<S+i.^h+i/2) - FM-i.<S, A6.-1/2)) , 

where as before q = (/i, hu^^, haxx, h<^zz): sjnd A6i_|_i/2 = ^i+i — bi. Thus we need to define the left 
and right numerical fluxes Fi{qi, qr, Ab), Fr{qi,qr, Ab), for all left and right values qi, qr, bi, br, 
with Ab — br~bi. We use the hydrostatic reconstruction method of P] (see also [II]), and define 

(5.52) h\ = {hi^{Ab)+)^, hl={hr-{-Ab)+)^, 

(5.53) qj = (j2f,hfu° i,hfaxx,i, h\(Tz2.,?j , ql = {hi, hlu° r, hlaxx,r, hla^.^ 
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with the notation .t+ = max(0, x). Note that we use the notation (1 instead of * in order to avoid 
confusions with the intermediate states of the Riemann solver of the previous sections. Then the 
numerical fuxes are defined by 

Fi{qi,q,., Ab) = Hq!,qt) + 0,.9^ -9^,0,0 

(5-54) > L . 

F,(g,,g,,Afo)=J-,(g?,qa)+ (^0,3^ -5^,0, Oj, 

where Ti and J> are the numerical fluxes (j5.40p of the problem without topography. 

Theorem 2. The scheme (|5.5ip with the numerical fluxes Fi, F^ defined by (|5.54p . (|5.52p . (|5.53p 

satisfies the following properties. 

(i) It is consistent with (j5.47p /or smooth solutions, 

(ii) It keeps the positivity of h, a^x, f^zz under the CFL condition AtA{qf,ql) < ^ min(Aa;;, Ax^) 
with A defined by (|5.44p . 

(Hi) It is conservative in the variable h, 

(iv) It satisfies a semi-discrete energy inequality associated to (|5.48p . 

(v) It is well-balanced, i.e. preserves the steady states at rest (|5.50p . 

Proof. We ommit the proof of the points (i) to (iii), which follow the proof of Proposition 4.14 in 

m- 

For the proof of (v), consider data g;, q^, bi, br at rest, i.e. satisfying ^ = = 0, hi+bi = hj.-\-br, 
o'xx,i = crxx,r = o-zz,i = o^zz.r = 1- Then from (j5.52p . (|5.53p we get qj = gj, the common value 

being q,. if A6 > 0, or qi if Ab < 0. We observe that then J'liqf.ql) = J'r{q!,ql) = F{q^) 
with F given in ([QS]) . and that indeed F{q'^) = {0, gh^^ /2, 0,0). The formulas yield Fi = 

(0,5/1^/2,0,0) = F{qi), Fr = (0,5/1^/2,0,0) = F{qr). If this is true at all interfaces, (|53T|) gives 
g"^^ ~ g.p, which proves the claim. 

Let us finally prove (iv). First, the scheme without topography satisfies the discrete energy 
inequality (|5.19p . According to [T2] section 2.2.2, it implies the semi-discrete energy inequality, 
characterized by 

ggs G{qr) + E'{qr){J^r{qi,qr) ^ F{qr)) < g{qi,qr), 

g{qi,qr)<G{qi)+E'{qi){Ti{qi,qr)-F{qi)), 

for all values of qi, qr, and where E' is the derivative of E with respect to q. Then, for the scheme 
with topography, the characterization of the semi-discrete energy inequality writes 

,5 5g^ Giqr, br) + E'{qr,br){Fr - F{qr)) < g{qi,qr, bl,br), 

g{qi,qr,bi,br)<G{qi,bi)+E'{qi,bi){Fi-F{qi)), 

where E and G are defined by (|5.49p . E' denotes the derivative of E with respect to q, and Q is 
an unknown consistent numerical entropy flux. Let us choose 

(5.57) g{qi,qr,bi,br) ^ Giqlqt) + J^'^qlqDgb^ 
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(5.61) 



where is the common /i— component of Ti and J>., and for some 5** that is defined below. Then, 
noticing that E'{q, b) = E'{q) + gb{l, 0, 0, 0), we can write the desired inequahties (|5.56p as 

(5 58) ^^'^'^^ + E'{qr){Fr - F{qr))+F\q\,ql)ghr < 5 (g?, (?».)+ J^^gf, 

^ ■ ^ g{qf,qt)+T''{q!,qt)gbi < G (qi) + E' (qi) {F, - F {qi)) + T'^ {qf , qt) gbi . 

But using (|5.55p evaluated at qf, and comparing the result with (|5.58p . we get the sufficient 

conditions 

(5.59) 

G(g,) + E'{qr){Fr - F{qr)) + T'^iqlqDgbr < G(g«) + i?'(g»)(J-,(gf , <?«) - F(g«)) + T''{qf,qt)gb\ 
G{qj) + E\ql){Ti{qf,qt) - F{ql)) + T\q\,ql)gb^ < G{qi) + E' {qi){Fi - F{qi)) + T\q\,ql)gbi. 

We compute now 

(5.60) E'{q) = (^-Ml! +gh-'-^ ln(a..a,,), ^(1 - 1/a..), |^(1 - 1/a,,) 

and writing 

F{q) = [hul h{ul)^ + g^ + ^/i(a,, - a,,,),ha,,,ulha,,u{^ , 
G{q) = + + + - ln('T,,a,,) - 2) + ^;i(a,, - a,,)^ u^, 

we deduce the identity 

(5.62) G{q)-E'{q)F{q)^~g^ul. 
Thus the inequality (j5.59p simplifies to 

^g^ul, + i?'(g,)F, + J^'^iqlqDghr < -g^u^, + E' {ql)Fr[q\,ql) + T\q\,ql)gbK 
(5-63) 2 2 2 

Now, using ([53^]) and the fact that E'{qr) - E'{ql) = {g{hr - /i«), 0,0,0), E'{qi) - E'{qj) = 
{g{hi — hf), 0,0,0), the desired inequalities (|5.63p rewrite 

5(/»r-^«+&.-6«)^"(gf,g«)<0, 
^ ^ ^ ff(/.,-;^r5» + 5,)-F^g?,g«)>0. 

We choose now 6" = max(6;, br), so that (|5.64p can be put in the form 

fr.r. {hr-ht-{-Ab) + )P\qf,ql)<0, 

^ ^ ^ {h,-hl-{AbU)T%lql)>0. 

Finally, taking into account (|5.52[) . we observe that if hi — (Afe)+ > then the second line of (j5.65|) 
is an identity, otherwise ~ and the the second inequality of (|5.65p holds because J^''(0, g|) < 
by the /i— nonnegativity of the numerical flux. The same argument is valid for the first inequality 
of (j5.65p . which concludes the proof. □ 
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Longitudinal stress convergence (50, 100, 200 and 400 points) 
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Figure 1 . Convergence of the discretized variables ha^x in Test case 1 



Remark 6. The maximum principle property on Sxx o-nd minimum principle property on s^z, that 
hold for the solver without topography, are not valid for the above solver with topography, even if 
it should hold at the continuous level. 



5.5. Numerical results. We now illustrate our model by numerical simulations performed with 
the scheme described above. Note that the model can be considered independently of its derivation 
and we explore numerical values beyond the physical regime of the Section [31 We denote by H{x) 
the Heaviside function with jump +1 at x = 0. For all numerical simulations, we chose Neumann 
conditions at boundary interfaces. 

Test case 1. It is a Riemann problem with initial condition {h, hu^, hoxx, hazz){t = 0) = 
(3 — 2H{x)){l, 0, 1, 1), without source term (6 = 0), that can be interpreted as a "dam" break on a 
wet floor, with polymeric fluid initially at rest everywhere. We first fix 77^ = A = 1 and study the 
convergence of our scheme with respect to the spatial discretization parameter for 50, 100, 200 and 
400 points and a constant CFL = 1/2. The results at final time T = .2 are shown in Fig.[TJ 
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Then, using 400 points and a constant CFL = 1/2, we let rjp vary as A = 1 is fixed. The 
results are shown in Fig. [5] and [3] On the one hand, the limit case 77p — J- coincides with the usual 
shallow- water model, since then the pressure assumes the same values as in a relaxation scheme for 
the Saint- Venant equations (independent of haxx and ha^z) while s^x^s^z become passive tracers 
and their evolution is only one-way coupled - in fact enslaved - to the autonomous dynamics 
of the Saint- Venant system of equations. On the other hand, the case j^p — >■ oo is some kind of 
"rigid limit". As expected from the formulae (|4.2ip for the eigenvalues of the Jacobian matrix, the 
left-going rarefaction wave and the right-going shock wave are all the faster as the viscosity rjp 
increases (we hardly see them at T = .2 for rjp = 10+'^ in Fig. [5] and [21). This is consistent with the 
physical notion of rigid limit (sound waves are faster in solids than in liquids) . On the contrary, 
the intermediate wave (a right-going contact discontinuity) is all the slower as rjp increases, and 
the jump of h across it is all the larger. This was not obvious to us at first. It could be explained 
for instance by the fact that ci -\- Cr (at the denominator of the formulae for u° ^ and the two 
intermediate states l/h^,l/h* in (|5.36p ) becomes very large when rjp increases. Notice that the 
jumps for axx and <7zz are directly related to that for h through (j5.38p and are consistently small 
when that for h is large. In any case, the materials on the left of the intermediate wave (where h is 
higher) is stretched in direction and compressed in direction , while it is stretched in direction 
Cz and compressed in direction on the right, which we can interpret as the manifestation of a 
close-to-equilibrium stability property (see Section |6]). 

We also let A vary as = 1 is fixed. The results (still at same given time T = .2) are shown in 
Fig. [Hand O One clearly sees here the competition between transport and diffusion of viscoelastic 
effects (the viscoelastic energy is stored in the new variables, that are transported but also diffused). 
Indeed, the viscoelastic energy is dissipated all the more rapidly as the relaxation time A is small 
(the viscous, "Low-Weissenberg" limit). And the waves are all the more smoothed as A is small 
(the source terms, which act as diffusive terms, are all the more important). On the other hand, 
the jump across the contact discontinuity is all the smaller for axx and (Tzz as A is small, but all the 
larger for h. This is coherent with a reasoning similar to the one above when rjp only was varied: 
a smaller relaxation time A implies slightly faster rarefaction and shock waves because of (|4.2ip . 
and a larger coefficient c/ -I- Cr in the formulae (|5.36p . 

Test case 2. It is a Riemann problem again without source term 6 = but with vacuum in 
the initial condition (h, hu'^, haxx, hazz)(t = 0) = (3 — 3if(a;))(l, 0, 1, 1), which can be interpreted 
as a "dam" break on a dry floor. The results in Fig. [6] and Fig. [7] at T — .5 show again that small 
A and large rjp imply a fast right-going rarefaction wave and a slow contact discontinuity, with a 
large jump for h and small jumps for axx,o'zz at contact discontinuity. On the contrary, large A 
and small rjp imply a slow right-going rarefaction wave and a fast contact discontinuity, with a 
small jump for h and large jumps for axx,<^zz at contact discontinuity. We nevertheless note that 
the "High-Weissenberg" limit A +oo (recall Remark [1]) is more difficult. In particular, although 
azz remains bounded, axx seems to become unbounded close to the wet/dry front. Recalling 
our Remark [5| (see also the comments below about the occurence of vacuum in Test 3), this is 
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Water height h for different r]p 




Figure 2. Variations of the variables with r/p in Test case 1 
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Longitudinal conformation a^x for different T]p 




Transversal conformation a^z for different ?/p 
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URE 3. Variations of the variables Oxxt^^zz with r/p in Test case 1 
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Water height h for different A 




not surprising: there hold minimum principles a^^ > {kh)~'^ ,azz > {khY in the High-Weissenberg 
limit A — > oo, ?/p/A = - without topography - where the relaxation source terms are neglected 

(except that here initially k — oo. but k becomes hopefully finite after some time). Hence a possible 
blow-up at the front, where — !• 0, when A is too large for the diffusive relaxation source terms 
to compensate for the transport effects. Note that in any case, the materials on the left of the 
contact-discontinuity wave (where h is non-zero) is stretched in direction Bx and compressed in 
direction e^,, except at the vacuum front where a small region with a closer to the equilibrium I 
is constantly seen (may be artificially due to our numerical treatment of the vacuum front). While 
no actual blow-up occurs, this will also be interpreted as a manifestation of a close-to-equilibrium 
stability property (see Section [6]). 

Test case 3. It is a benchmark introduced in [33] to test the treatment of topography, see 
also d^. For x G (0,25), we compute until T = .25 the evolution from an initial condition 
(/i, hul, ha^x, hazz){t 0, x) = ((10 - 6)+, -350 + 700H{x - 50/3), (10 - 6)+, (10 - b)+) over a 
topography b{x) ~ H{x — 25/3) — H{x — 25/2). Two rarefaction waves propagate on the left and 
right sides of the initial velocity singularity at x = 50/3 so that a vacuum is created in between 
(in the usual Saint- Venant case). In addition, a couple of rarefaction/shock waves is created at 
each singular point x = 25/3, 25/2 of the topography, but have much smaller amplitudes than the 
rarefaction waves at x = 50/3. The results in Fig. [5] and [5] are obtained for various rjp,\ at a 
constant rjp/X = lO""'. This particular choice was made because then, at the final time, the system 
is sufficiently close to the Saint- Venant limit ?/p/A so that the pressure is hardly modified 
compared with the usual Saint- Venant case. 

Compared with the usual Saint- Venant case, the double rarefaction wave centered at a: = 50/3 
cannot create vacuum but at the single point x = 50/3 where the initial velocity is singular. This 
can be explained as follows. Assuming that the source terms in the stress equations do not influence 
much the bounds on axx,(^zzi in agreement with our Remark 21 the maximum (resp. minimum) 
principle holds for Sxx (resp. Szz), and there exists a constant k (depending only on the initial 
conditions since initially ft. > 0) such that Gxx > {kh)~'^ , azz > {kh)^'^ . But according to the 
energy bound, one has that (rip/X) J haxxdx remains bounded. We deduce that {r]p/X) J h~^dx 
remains bounded, and therefore h cannot tend to on a whole interval, but can vanish on a single 
point. We have then cFxx ^ +oo at the singular point, here x ~ 50/3. 

Moreover, another vacuum can be created a.t x — 25/2 (of course still at a single-point because 
of the previous reasoning). Of course, this new phenomena could just be a numerical artifact, for 
instance due to the computation of source terms by the hydrostatic reconstruction method. In 
any case, assuming it is part of the richer phenomenology of our new model compared with the 
usual Saint- Venant model (see also Section [5] for possible physical interpretations), we observe that 
the existence of that phenomena depends on yyp, A. It happens only for rip ~ 10^**, A = 10+" and 
r\p = 10~^,A = 10+^ in our numerical experiments, not for r]p = 10~^,A = 10^^. But it occurs 
for larger A at i]p = 10~^ (when the stress relaxation terms are less important). Indeed, this 
phenomena seems triggered mainly by high values of A (the "High-Weissenberg limit", which by the 
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Water height h for different r/p , A 




Water velocity w° for different rjp, X 
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Figure 6. Variations of the variables h, u'^ with r/p, A in Test case 2 
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Longitudinal conformation Uxx for different ?7p,A 
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Longitudinal conformation azz for different ?7p,A 
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Figure 7. Variations of the variables cr^x, cr^^ with rjp, A in Test case 2 



way requires very small time steps because of the CFL constraint), and still holds at high values of 
rjp (a "rigid limit" that seems to lead to some kind of a "break-up" of the jet here) , even at very large 
r]p when no vacuum occurs in between the rarefaction waves. Notice that the latter phenomenon 
also induces an additional sign change for in between the two vacuum points, whose location 
seems to depend on A but not on r/p, (another indication that the role of viscoelastic dissipation 
is dominant here). And compared with the usual Saint- Venant case, velocities also assume much 
greater value (on the left in particular). 

Test case 4. In our last test case, we woud like to assess the treament of another type of 
topography source terms, with creation of dry/wet fronts, by the hydrostatic reconstruction. A 
usual test case is Thacker's [46] e.g., which has analytical solutions in the usual Saint- Venant case. 
But we could not capture interesting phenomena with Thacker's testcase, in particular because 
the CFL constraint requires the time step to go to very quickly (on short time ranges) due to 
the creation of dry fronts where /i — >■ and (j^x +oo (with A not too small) . Note that this does 
not necessarily mean that this problem does not have global solutions with finite-energy. A time- 
implicit scheme (probably hard to build) might be able to compute finite-energy approximations 
with a non- vanishing time-step. 

Here, we consider the test case proposed by Synolakis |1S] to model the runup of solitary 
waves. This could be used until interesting final times T = 32.5 after the incidental wave has 
refiected against the shore and created a dry front, see also [38|. We use {h,hu'^,haxx,hazz){t = 
0,x) = ((1. + hQ{x) — b{x))^ (1, Y^/io(a;), 1, 1) as initial condition over a topography b{x) ~ {{x — 
40.)/19.85)+, X G (0,100). Thepertubation/io(x) = a{cosh{y/j5^{x~acosh{,/T/^)/{.75a))))~'^ 
models a solitary wave as a function of the parameter with a = .019/.1 according to Synolakis 
semi-analytical theory. 

The results in Fig. [TUl and [TT] show that it is essentially the variations of r/p/X that infiuence 
the water height and velocity among all possible variations of r/p, A. And although the first effect 
of the variations of rjp/X is on the waves celerity, there is no direct match between variations in 
r]p/X and a time shift as shown in Fig. 1101 On the contrary, the variables axx^f^zz depend more 
on A alone (recall the importance of relaxation source terms, especially in the case where /i — > 0), 
at least for such small values of rjp/X as those tested here (sufficiently close to the Saint- Venant 
regime for the time step not to vanish, even at high values of A). The smaller A is, the stronger the 
dissipation is and thus balances the large stress values that were induced close to the dry front. 

6. Conclusion 

We have proposed a new reduced model for the motion of thin layers of viscoelastic fluids 
(shallow viscoelastic flows) that are described by the upper-convected Maxwell model and driven 
by the gravity, under a free surface and above a given topography with small slope (like in the 
standard Saint Venant model for shallow water). More precisely, we have shown formally that 
for given boundary conditions and under scaling assumptions (Hl-5), the solution to the incom- 
pressible Euler-UCM system of equations can be approximated by the solutions to the reduced 
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Figure 8. Variations of the variable h + b,u'^ with r/p, A in Test case 3. We use 
different labels for the positive (©) and negative (g) part of the velocity. 
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Longitudinal conformation (Tj^^ for different r/p, A 
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Longitudinal conformation cr^^ for different rjp, A 
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Figure 9. Variations of the variable axx,crzz with ?/p, A in Test case 3. 
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Water height h for different 7]p , X 
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Water velocity m° for different r/p, A 
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Figure 10. Variations of the variable h + b,u^ with r]p,X in Test case 4 
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Longitudinal conformation cr^:!: for different rjp , A 
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Longitudinal conformation cr^^ for different rjp, A 
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Figure 11. Variations of the variable axx,crzz with rip,X in Test case 4 
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model (j3.22p in some asymptotic regime. We hope that this asymptotic regime in particular is 
physically meaningful, that and our new model makes sense, possibly beyond the previous asymp- 
totic regime. (That is why we have studied it mathematically and explored it numerically without 
constraining ourself to a particular regime, as it is usual in such cases.) 

Observe that in the end we have obtained a flow model whose dynamics is function of the first 
normal stress difference only, while the shear part of the stress is negligible and computed as an 
output of the flow evolution. More specifically, the boundary conditions (|2.10| - f^?5|) and the flat 
velocity profile (consequence of the assumed motion by slices) require compatibility conditions on 
the bulk behaviour of t^z inside a thin layer. Before looking in future works for other asymptotic 
regimes, possibly compatible (under different assumptions) with more general kinematics, we would 
like to conclude here with a better insight of the physical implications of our reduced model. 

6.1. Physical interpretation from the macroscopic mechanical viewpoint. We note that 
the main differences between our model for shallow (Maxwell) viscoelastic flows and the standard 
Saint Venant model for shallow water is i) a new hydrostatic pressure p.lOp . which is function of 
the (viscoelastic) internal stresses in addition to the water level h, hence ii) a new hydrodynamic 
force in the momentum balance (in addition to the external gravity force), which is proportional 
to the normal stress difference t^x — '''zz^ and iii) variable internal stresses t^x and r^^, which have 
their own dynamics corresponding to a viscoelastic mechanical behaviour (with a finite relaxation 
time A = 0(1) ; such that one recovers the standard viscous mechanical behaviour only in the limit 
A — > 0). Moreover, in the asymptotic regime where our non- Newtonian model was derived, with a 
small viscosity parameter r\p = 0(e), the strain and stress tensors have the scaling 



One essential rheological feature of our reduced model is thus the ratio e between the shear and 
elongational components of the stress tensor r and of the strain tensor Vii. The fact that our 
model should mainly describe extensional flows, with small shear (of the same small order as the 
elongational viscosity), seems to be a strong limitation to the applicability of our model in real 
situations. Of course, one is likely to need another reduced model (in other asymptotic regimes) 
to describe flows that are not essentially elongational. 

Note yet that there are situations where physicists arrive at similar conclusions [221 [23] and 
obtain a very similar one-dimensional model with purely elongational stresses for the description of 
free axisymmetric jets. By the way, a description of free axisymmetric jets is also well achieved by 
our model since the pure slip boundary conditions (|2.7p - (|2.8p is equivalent to assuming a cylindrical 
symmetry around the symmetry line of the jet, and surface tension effects (neglected in our model) 
can be included using standard modifications of our no-tension boundary condition (|2.10p . 

Moreover, it seems possible to still include non-negligible shear effects in our model through a 
parabolic correction of the vertical profile like in 137). as well as surface tension and friction 
effects of order one at the boundaries. 



(6.1) 
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6.2. Physical interpretation at the microscopic molecular level. A microscopic interpre- 
tation of our asymptotic regime can also be achieved using a molecular model of the elastic effects 
(that is, a model at the molecular level from which the UCM is a coarse-grained version at the 
macroscopic mechanical level). Following [TU], a typical molecular model that accounts for the 
elasticity of a fluid invokes the transport of elastically deformable Brownian particles diluted in 
the fluid (which can often be thought of as large massive molecules like polymers). The simplest 
model of this kind couples, locally in the physical space, a kinetic theory for "dumbbells" (two 
point-masses connected by an elastic force idealized as a "spring") with the strain of the fluid. 

Let us denote Xt{x) the connector vector between the two point- masses of a dumbbell modelling 
a polymer molecule at position x and time t in the fluid. The collection of vector stochastic 
processes {Xt{x))f^(^Q^^^-^ parametrized by a; G 2?4 is solution to overdamped Langevin equations 



(6.2) dXf + {u ■ V)Xtdt = [{Vu)Xt - |F(Xt)^ dt + 2^ ^dB* 

for a given field {Qt{x))t^{o^+ao) of standard Brownian motions (in Ito sense) where C is a friction 
parameter, ks the Boltzmann constant and T the absolute thermodynamical temperature. The 
UCM equations can be exactly recovered with the specific choice F(X) = HX. Indeed, the 
extra-stress t and the conformation tensor cr are given by Kramers relation 

(6.3) T = ^(TOcr-/) with a{t,x) = ^K\Xt®¥{Xt)] = I \x ® x]'tjj{t,x,iX)dX 

2A Hi 1 J I 1 

where X = X/^ is an adimensional version of X, to = is a ratio between the elastic potential 
energy and the heat of the Brownian bath, and rjp = 2XnkBT is the molecular interpretation of 
the polymer viscosity, with n the number density of polymer chains by unit volume (assumed 
constant as usual for dilute polymer solutions) and A a characteristic time for dumbbells. One can 
always choose i such that m = 1. Then, choosing A = ^ as a relaxation time, Ito formula allows 
one to exactly recover the UCM system of equations (|2.3p when the solvent is assumed inviscid 
with a velocity field u{t, x) solution to the Euler equations, on noting that the probability density 
ip{t, x,£X) satisfies the following Fokker-Planck equation on the unbounded domain X e M^: 

(6.4) ^+^,. vv,' = -div^ (^[(V«)X-^X]^) +^A^V. 

We note that in [1^ a reduced model for shallow viscoelastic flows quite similar to ours has 
already been derived starting from a coupled micro- macro system like (|6.4H6.3H2.2p . rather than 
starting from a coarse-grained system at the macroscopic level like the UCM model. The difference 
between the Hookean micro-macro system above (equivalent in some sense to the UCM model) 
and the micro-macro system used in [JD] is the spring force: it corresponds to FENE dumbbells 
F(Xt) = Xt/(1 — |Xtp/6) in [3D]. The FENE force is more physical because it accounts for a 
finite extension |Xi| < b, but contrary to the Hookean dumbbells, it does not have an exact coarse- 
grained macroscopic equivalent like the UCM model. Yet, if we follow the same procedure as in [3D] 
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but for Hookean dumbbells, we can hope to derive a reduced micro-macro model whose coarse- 
grained version is comparable to our new reduced UCM model. Moreover, if the scaling regimes 
are the same as in |3D], then our model should also compare to that in |3D]; for an inviscid solvent, 
in the infinite extensibility limit b ^ oo (where one formally recovers the Hookean dumbbells from 
FENE dumbbells). Now, observe that the scaling of our new model implies (j6.ip Vit = 7o + 0{e) 
where 7q = 0(1) is a traceless diagonal matrix with entries dxu'^, —dxu'^- Then (|6.4p rewrites 

where M{t,x, X) is a weight function proportional to the Maxwellian e~-^ i^^ta^^)^ ^ xhe ap- 
proximation (|6.5p of (j6.4[) is consistent with our new reduced model provided it yields a consistent 
approximation for the stress in (j6.3p : that is, it suffices to show axxT^Tzz = 0{1) and axz = 0{e) 
as e — >■ 0. To this aim. let us define an order-one approximation xp'^ = ijj + 0(e) solution to 

(6.6, f^ + „..v*« = i;d.v.(MV.(|; 

The point is to estimate the terms 

(6.7) r° = |^(^t"-7), ct" = j [X<»X]^°{X)dX . 

This is not an easy task because of the coupling between and it°. Yet it seems reasonable 
to assume that V'^ remains close to the equilibrium solution M/ J M for all times (indeed, the 
Hookean force is derived from an a-convex potential [2]); and in particular the Maxwellian ilP oc 
^-{Ax^+Bz^-+2Cxz) j^j^g ^j^g scaling A = 0{l),B = 0(1),C = 0(e), which implies that a (and thus 
r) is diagonal at first order. One then obtains a reduced kinetic model which can be exactly 
coarse-grained into our new reduced UCM model with Ito formula. 

A macroscopic consequence of the microscopic assumption above is that the reduced model is 
well-adapted for elongational flows, which is consistent with our macroscopic intepretation of the 
model. Indeed, everywhere in the macroscopic physical space, one can only expect a balance of 
internal elastic energy due to stretching or compressing strains in the directions Bx and , which is 
the case in elongational flows. (There is not a high probability of permanently sheared dumbbells.) 
Moreover, if ifp is actually close to the equilibrium M (the particular case A k, 2XdxU^ — 1,_B w 
—2XdxU^ — 1,C sa of our assumption), then, at first-order, the dumbbells are quite uniformly 
oriented but stretched in one canonical direction - Bx or - and necessarily compressed in the 
orthogonal one (the level-sets of the distribution function are ellipsoidal with principal axes Sx 
and Bz at first order) . This was indeed observed in those numerical experiments where no blow-up 
phenomenon seemed to occur. 

The microscopic view is in turn a plausible physical explanation at the molecular level of some 
macroscopic observations. Recall indeed that one-dimensional simple models similar to our model 
have already been derived in the past to model axisymmetric free jets of elastic liquids [551 123 
with a view to explaining the die swell at the end of an extrusion pipe. Now, a miscroscopic 
interpretation of the die swell is: the elastic energy stored before the die is released after the die. 
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The dumbbells, mainly compressed in the radial direction before the die, stretch just after 
the die. This may be responsible for an increase of the jet radius (the free-surface of the jet flow 
equilibrates with the atmospheric pressure) after a characteristic relaxation time linked to A, hence 
the so-called delayed die swell. 

Finally, we would like to comment on the results obtained in |3D] with FENE dumbbells. The 
main differences with our reduced model (which has the micro-macro interpretation explicited 
above) are: (i) the relaxation time in [ID] is assumed small A = 0(e), because then it is possible 
to compute approximate solutions to the Fokker-Planck equation following the Chapman-Enskog 
procedure of |21| ; and (ii) the polymer distribution is mainly radial {ip^ is only function of \X.\), 
because the authors claim that this suffices to next imply cixz = 0{e) and, as a consequence, a 
flat proflle for the horizontal velocity like in our model. Then, the scaling regimes are not the 
same, and the radial assumption is too strong to allow one to recover our ellipsoidal probability 
distribution. So we cannot directly compare our results though they have a similar flavour. 



6.3. Open questions and perspectives. First, regarding the interpretation of our model, one 
might ask whether the present scaling corresponds to a physical situation actually observed for 
elastic fluids in nature. In particular, the main questionable assumption is of course the pure-slip 
and no-friction boundary conditions (j2.7fl2.8p at the bottom (already unrealistic for Newtonian 
flows, maybe even more unrealistic for non-Newtonian ones). Second, future works on this topic 
might consider the following directions: 

• derive thin-layer reduced models with other equations modelling non-Newtonian flows, 
which are believed to better model the rheological properties of real materials (constitutive 
models like Giesekus, PTT, FENE-P, or other molecular models than the FENE dumbbell 
model used in [111 HO]), and in two-dimensional settings (see [151 ISZ] for the standard 
shallow water model); 

• derive a reduced model closer to real physical situations, possibly in different regimes, or 
for instance by using a ^i-dependent velocity proflle u.^ (possibly a multi-layer model) and 
different boundary conditions than (|2.10p and (|2.8p (with surface tension and friction at 
the bottom), which may lead to flnd physical regimes where t^z is not negligible; 

• give a rigorous mathematical meaning and enhance numerical simulations (well-balanced 
second-order reconstructions) for non-standard systems of equations like the new one pre- 
sented here. 

We note that multi-layer models are also a path to the modelling of some important physical 
situations, like a thin layer of polymeric fluids on water to forecast the efficiency of oil slick 
protection plans. 
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Appendix A. Convexity of the energy 



In order to check the convexity of E in (|4.14p with respect to general variables, we use a 
Lagrange transformation, see for example Lemma 1.4 in [T^]. Thus E is a convex function of 



h, hu'^, hw 



-1/2^ 



, hi; 



' 1/2^ 



for given smooth invertible functions tn, if and only if iJ/ft, is a convex function of the Lagrangian 
variables 



' -1/2^ 



1/2^ 



Let us denote by V^, i = 1, . . . , 4 the entries of the vector then the Lagrangian energy writes 



E 



1 9 



9 1 



9h 



In 



- 2 



h 2' 2^1 4A 1^^(1/3)2 -y^^y^Y, , 

Litroduce now the notation 

n{V^) = 2\nm{V:i), C(^4) = -21n<j(F4). 

Clearly we only need to look at the convexity with respect to Vij, V4), and the Hessian matrix 
T-L oi E /h with respect to these variables (at fixed h) is given by 



■Tip 



^^ + 2e-- + ^ 



2^ 



-2Vie-^^Vt' 








^(C'2-C")+C" 

where prime denotes the derivative with respect to the involved Vi . Since Vi can take any positive 
value at fixed V3 or V4, the positivity of the diagonal terms give the necessary conditions 



< o"(y3) < ^'{v^)\ < C'iVi) < CiVi)^- 

Then, writing the positivity of the determinant of the 2x2 upper left submatrix of "H, and looking 
at the dominant term when Vi — > 00 yields the necessary condition 



2e 



-20^ri'2 



Obviously there is no function ^{Vs) satisfying these conditions, and E is never convex with respect 
to the considered variables. 

On the contrary, if we choose the physically natural, but non-conservative, variables q — 
{h, hu'^, huxxi hazz), then using the Lagrangian variables 

1 



W 



, W^, (Jxxi ^zz I : 



one can write 



E 
'h 



+ i^+gb+^ {axx + (Tzz - In {(Txx<^zz) - 2) 
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which is obviously convex with respect to W (at fixed b). We conclude that E is convex with 
respect to q. 
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